//----------------------------------------------------------------------
//	File:		kd_split.cc
//	Programmer:	Sunil Arya and David Mount
//	Last modified:	03/04/98 (Release 0.1)
//	Description:	Methods for splitting kd-trees
//----------------------------------------------------------------------
// Copyright (c) 1997-1998 University of Maryland and Sunil Arya and David
// Mount.  All Rights Reserved.
// 
// This software and related documentation is part of the 
// Approximate Nearest Neighbor Library (ANN).
// 
// Permission to use, copy, and distribute this software and its 
// documentation is hereby granted free of charge, provided that 
// (1) it is not a component of a commercial product, and 
// (2) this notice appears in all copies of the software and
//     related documentation. 
// 
// The University of Maryland (U.M.) and the authors make no representations
// about the suitability or fitness of this software for any purpose.  It is
// provided "as is" without express or implied warranty.
//----------------------------------------------------------------------

#include "ANN.h"
#include "ANNx.h"			// all ANN includes
#include "ANNperf.h"		// performance evaluation

#include "kd_tree.h"			// kd-tree definitions
#include "kd_util.h"			// kd-tree utilities
#include "kd_split.h"			// splitting functions

//----------------------------------------------------------------------
//  Constants
//----------------------------------------------------------------------

const double ERR = 0.001;		// a small value

//----------------------------------------------------------------------
//  kd_split - Bentley's standard splitting routine for kd-trees
//	Find the dimension of the greatest spread, and split
//	just before the median point along this dimension.
//----------------------------------------------------------------------

void kd_split(
    ANNpointArray	pa,		// point array (permuted on return)
    ANNidxArray		pidx,		// point indices
    const ANNorthRect	&bnds,		// bounding rectangle for cell
    int			n,		// number of points
    int			dim,		// dimension of space
    int			&cut_dim,	// cutting dimension (returned)
    ANNcoord		&cut_val,	// cutting value (returned)
    int			&n_lo)		// num of points on low side (returned)
{
					// find dimension of maximum spread
    cut_dim = annMaxSpread(pa, pidx, n, dim);
    n_lo = n/2;				// median rank
					// split about median
    annMedianSplit(pa, pidx, n, cut_dim, cut_val, n_lo);
}

//----------------------------------------------------------------------
//  midpt_split - midpoint splitting rule for box-decomposition trees
//
//	This is the simplest splitting rule that guarantees boxes
//	of bounded aspect ratio.  It simply cuts the box with the
//	longest side through its midpoint.  If there are ties, it
//	selects the dimension with the maximum point spread.
//
//	WARNING: This routine (while simple) doesn't seem to work
//	well in practice in high dimensions, because it tends to
//	generate a large number of trivial and/or unbalanced splits.
//	Either kd_split(), sl_midpt_split(), or fair_split() are
//	recommended, instead.
//----------------------------------------------------------------------

void midpt_split(
    ANNpointArray	pa,		// point array
    ANNidxArray		pidx,		// point indices (permuted on return)
    const ANNorthRect	&bnds,		// bounding rectangle for cell
    int			n,		// number of points
    int			dim,		// dimension of space
    int			&cut_dim,	// cutting dimension (returned)
    ANNcoord		&cut_val,	// cutting value (returned)
    int			&n_lo)		// num of points on low side (returned)
{
    int d;

    ANNcoord max_length = bnds.hi[0] - bnds.lo[0];
    for (d = 1; d < dim; d++) {		// find length of longest box side
	ANNcoord length = bnds.hi[d] - bnds.lo[d];
	if (length  > max_length) {
	    max_length = length;
	}
    }
    ANNcoord max_spread = -1;		// find long side with most spread
    for (d = 0; d < dim; d++) {
					// is it among longest?
	if (double(bnds.hi[d] - bnds.lo[d]) >= (1-ERR)*max_length) {
					// compute its spread
	    ANNcoord spr = annSpread(pa, pidx, n, d);
	    if (spr > max_spread) {	// is it max so far?
		max_spread = spr;
		cut_dim = d;
	    }
	}
    }
					// split along cut_dim at midpoint
    cut_val = (bnds.lo[cut_dim] + bnds.hi[cut_dim]) / 2;
					// permute points accordingly
    int br1, br2;
    annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
    //------------------------------------------------------------------
    //	On return:	pa[0..br1-1] < cut_val
    //			pa[br1..br2-1] == cut_val
    //			pa[br2..n-1] > cut_val
    //
    //	We can set n_lo to any value in the range [br1..br2].
    //	We choose split so that points are most evenly divided.
    //------------------------------------------------------------------
    if (br1 > n/2) n_lo = br1;
    else if (br2 < n/2) n_lo = br2;
    else n_lo = n/2;
}

//----------------------------------------------------------------------
//  sl_midpt_split - sliding midpoint splitting rule
//
//	This is a modification of midpt_split, which has the nonsensical
//	name "sliding midpoint".  The idea is that we try to use the
//	midpoint rule, by bisecting the longest side.  If there are
//	ties, the dimension with the maximum spread is selected.  If,
//	however, the midpoint split produces a trivial split (no points
//	on one side of the splitting plane) then we slide the splitting
//	(maintaining its orientation) until it produces a nontrivial
//	split.  For example, if the splitting plane is along the x-axis,
//	and all the data points have x-coordinate less than the x-bisector,
//	then the split is taken along the maximum x-coordinate of the
//	data points.
//
//	Intuitively, this rule cannot generate trivial splits, and
//	hence avoids midpt_split's tendency to produce trees with
//	a very large number of nodes.
//
//----------------------------------------------------------------------

void sl_midpt_split(
    ANNpointArray	pa,		// point array
    ANNidxArray		pidx,		// point indices (permuted on return)
    const ANNorthRect	&bnds,		// bounding rectangle for cell
    int			n,		// number of points
    int			dim,		// dimension of space
    int			&cut_dim,	// cutting dimension (returned)
    ANNcoord		&cut_val,	// cutting value (returned)
    int			&n_lo)		// num of points on low side (returned)
{
    int d;

    ANNcoord max_length = bnds.hi[0] - bnds.lo[0];
    for (d = 1; d < dim; d++) {		// find length of longest box side
	ANNcoord length = bnds.hi[d] - bnds.lo[d];
	if (length  > max_length) {
	    max_length = length;
	}
    }
    ANNcoord max_spread = -1;		// find long side with most spread
    for (d = 0; d < dim; d++) {
					// is it among longest?
	if ((bnds.hi[d] - bnds.lo[d]) >= (1-ERR)*max_length) {
					// compute its spread
	    ANNcoord spr = annSpread(pa, pidx, n, d);
	    if (spr > max_spread) {	// is it max so far?
		max_spread = spr;
		cut_dim = d;
	    }
	}
    }
					// ideal split at midpoint
    ANNcoord ideal_cut_val = (bnds.lo[cut_dim] + bnds.hi[cut_dim])/2;

    ANNcoord min, max;
    annMinMax(pa, pidx, n, cut_dim, min, max);	// find min/max coordinates

    if (ideal_cut_val < min)		// slide to min or max as needed
	cut_val = min;
    else if (ideal_cut_val > max)
	cut_val = max;
    else
	cut_val = ideal_cut_val;

					// permute points accordingly
    int br1, br2;
    annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
    //------------------------------------------------------------------
    //	On return:	pa[0..br1-1] < cut_val
    //			pa[br1..br2-1] == cut_val
    //			pa[br2..n-1] > cut_val
    //
    //	We can set n_lo to any value in the range [br1..br2] to satisfy
    //	the exit conditions of the procedure.
    //
    //	if ideal_cut_val < min (implying br2 >= 1),
    //		then we select n_lo = 1 (so there is one point on left) and
    //  if ideal_cut_val > max (implying br1 <= n-1),
    //		then we select n_lo = n-1 (so there is one point on right).
    //	Otherwise, we select n_lo as close to n/2 as possible within
    //		[br1..br2].
    //------------------------------------------------------------------
    if (ideal_cut_val < min) n_lo = 1;
    else if (ideal_cut_val > max) n_lo = n-1;
    else if (br1 > n/2) n_lo = br1;
    else if (br2 < n/2) n_lo = br2;
    else n_lo = n/2;
}

//----------------------------------------------------------------------
//  fair_split - fair-split splitting rule
//
//	This is a compromise between the kd-tree splitting rule (which
//	always splits data points at their median) and the midpoint
//	splitting rule (which always splits a box through its center.
//	The goal of this procedure is to achieve both nicely balanced
//	splits, and boxes of bounded aspect ratio.
//
//	A constant FS_ASPECT_RATIO is defined.  Given a box, those sides
//	which can be split so that the ratio of the longest to shortest
//	side does not exceed ASPECT_RATIO are identified.  Among these
//	sides, we select the one in which the points have the largest
//	spread.  We then split the points in a manner which most evenly
//	distributes the points on either side of the splitting plane,
//	subject to maintaining the bound on the ratio of long to short
//	sides.  To determine that the aspect ratio will be preserved,
//	we determine the longest side (other than this side), and
//	determine how narrowly we can cut this side, without causing the
//	aspect ratio bound to be exceeded (small_piece).
//
//	This procedure is more robust than either kd_split or midpt_split,
//	but is more complicated as well.  When point distribution is
//	extremely skewed, this degenerates to midpt_split (actually
//	1/3 point split), and when the points are most evenly distributed,
//	this degenerates to kd-split.
//----------------------------------------------------------------------

const double FS_ASPECT_RATIO = 3.0;	// maximum allowed aspect ratio
					// must be >= 2.

void fair_split(
    ANNpointArray	pa,		// point array
    ANNidxArray		pidx,		// point indices (permuted on return)
    const ANNorthRect	&bnds,		// bounding rectangle for cell
    int			n,		// number of points
    int			dim,		// dimension of space
    int			&cut_dim,	// cutting dimension (returned)
    ANNcoord		&cut_val,	// cutting value (returned)
    int			&n_lo)		// num of points on low side (returned)
{
    int d;
    ANNcoord max_length = bnds.hi[0] - bnds.lo[0];
    cut_dim = 0;
    for (d = 1; d < dim; d++) {		// find length of longest box side
	ANNcoord length = bnds.hi[d] - bnds.lo[d];
	if (length  > max_length) {
	    max_length = length;
	    cut_dim = d;
	}
    }

    ANNcoord max_spread = 0;		// find legal cut with max spread
    cut_dim = 0;
    for (d = 0; d < dim; d++) {
	ANNcoord length = bnds.hi[d] - bnds.lo[d];
					// is this side midpoint splitable
					// without violating aspect ratio?
	if (((double) max_length)*2.0/((double) length) <= FS_ASPECT_RATIO) {
					// compute spread along this dim
	    ANNcoord spr = annSpread(pa, pidx, n, d);
	    if (spr > max_spread) {	// best spread so far
		max_spread = spr;
		cut_dim = d;		// this is dimension to cut
	    }
	}
    }

    max_length = 0;			// find longest side other than cut_dim
    for (d = 0; d < dim; d++) {
	ANNcoord length = bnds.hi[d] - bnds.lo[d];
	if (d != cut_dim && length > max_length)
	    max_length = length;
    }
					// consider most extreme splits
    ANNcoord small_piece = max_length / FS_ASPECT_RATIO;
    ANNcoord lo_cut = bnds.lo[cut_dim] + small_piece;// lowest legal cut
    ANNcoord hi_cut = bnds.hi[cut_dim] - small_piece;// highest legal cut

    int br1, br2;
					// is median below lo_cut ?
    if (annSplitBalance(pa, pidx, n, cut_dim, lo_cut) >= 0) {
	cut_val = lo_cut;		// cut at lo_cut
	annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
	n_lo = br1;
    }
					// is median above hi_cut?
    else if (annSplitBalance(pa, pidx, n, cut_dim, hi_cut) <= 0) {
	cut_val = hi_cut;		// cut at hi_cut
	annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
	n_lo = br2;
    }
    else {				// median cut preserves asp ratio
	n_lo = n/2;			// split about median
	annMedianSplit(pa, pidx, n, cut_dim, cut_val, n_lo);
    }
}

//----------------------------------------------------------------------
//  sl_fair_split - sliding fair split splitting rule
//
//	Sliding fair split is a splitting rule that combines the
//	strengths of both fair split with sliding midpoint split.
//	Fair split tends to produce balanced splits when the points
//	are roughly uniformly distributed, but it can produce many
//	trivial splits when points are highly clustered.  Sliding
//	midpoint never produces trivial splits, and shrinks boxes
//	nicely if points are highly clustered, but it may produce
//	rather unbalanced splits when points are unclustered but not
//	quite uniform.
//
//	Sliding fair split is based on the theory that there are two
//	types of splits that are "good": balanced splits that produce
//	fat boxes, and unbalanced splits provided the cell with fewer
//	points is fat.
//
//	This splitting rule operates by first computing the longest
//	side of the current bounding box.  Then it asks which sides
//	could be split (at the midpoint) and still satisfy the aspect
//	ratio bound with respect to this side.  Among these, it selects
//	the side with the largest spread (as fair split would).  It
//	then considers the most extreme cuts that would be allowed by
//	the aspect ratio bound.  This is done by dividing the longest
//	side of the box by the aspect ratio bound.  If the median cut
//	lies between these extreme cuts, then we use the median cut.
//	If not, then consider the extreme cut that is closer to the
//	median.  If all the points lie to one side of this cut, then
//	we slide the cut until it hits the first point.  This may
//	violate the aspect ratio bound, but will never generate empty
//	cells.  However the sibling of every such skinny cell is fat,
//	and hence packing arguments still apply.
//
//----------------------------------------------------------------------

const double SFS_ASPECT_RATIO = 3.0;	// maximum allowed aspect ratio
					// must be >= 2.

void sl_fair_split(
    ANNpointArray	pa,		// point array
    ANNidxArray		pidx,		// point indices (permuted on return)
    const ANNorthRect	&bnds,		// bounding rectangle for cell
    int			n,		// number of points
    int			dim,		// dimension of space
    int			&cut_dim,	// cutting dimension (returned)
    ANNcoord		&cut_val,	// cutting value (returned)
    int			&n_lo)		// num of points on low side (returned)
{
    int d;
    ANNcoord min, max;			// min/max coordinates
    int br1, br2;			// split break points

    ANNcoord max_length = bnds.hi[0] - bnds.lo[0];
    cut_dim = 0;
    for (d = 1; d < dim; d++) {		// find length of longest box side
	ANNcoord length = bnds.hi[d] - bnds.lo[d];
	if (length  > max_length) {
	    max_length = length;
	    cut_dim = d;
	}
    }

    ANNcoord max_spread = 0;		// find legal cut with max spread
    cut_dim = 0;
    for (d = 0; d < dim; d++) {
	ANNcoord length = bnds.hi[d] - bnds.lo[d];
					// is this side midpoint splitable
					// without violating aspect ratio?
	if (((double) max_length)*2.0/((double) length) <= SFS_ASPECT_RATIO) {
					// compute spread along this dim
	    ANNcoord spr = annSpread(pa, pidx, n, d);
	    if (spr > max_spread) {	// best spread so far
		max_spread = spr;
		cut_dim = d;		// this is dimension to cut
	    }
	}
    }

    max_length = 0;			// find longest side other than cut_dim
    for (d = 0; d < dim; d++) {
	ANNcoord length = bnds.hi[d] - bnds.lo[d];
	if (d != cut_dim && length > max_length)
	    max_length = length;
    }
					// consider most extreme splits
    ANNcoord small_piece = max_length / SFS_ASPECT_RATIO;
    ANNcoord lo_cut = bnds.lo[cut_dim] + small_piece;// lowest legal cut
    ANNcoord hi_cut = bnds.hi[cut_dim] - small_piece;// highest legal cut
					// find min and max along cut_dim
    annMinMax(pa, pidx, n, cut_dim, min, max);
					// is median below lo_cut?
    if (annSplitBalance(pa, pidx, n, cut_dim, lo_cut) >= 0) {
	if (max > lo_cut) {		// are any points above lo_cut?
	    cut_val = lo_cut;		// cut at lo_cut
	    annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
	    n_lo = br1;			// balance if there are ties
	}
	else {				// all points below lo_cut
	    cut_val = max;		// cut at max value
	    annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
	    n_lo = n-1;
	}
    }
					// is median above hi_cut?
    else if (annSplitBalance(pa, pidx, n, cut_dim, hi_cut) <= 0) {
	if (min < hi_cut) {		// are any points below hi_cut?
	    cut_val = hi_cut;		// cut at hi_cut
	    annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
	    n_lo = br2;			// balance if there are ties
	}
	else {				// all points above hi_cut
	    cut_val = min;		// cut at min value
	    annPlaneSplit(pa, pidx, n, cut_dim, cut_val, br1, br2);
	    n_lo = 1;
	}
    }
    else {				// median cut is good enough
	n_lo = n/2;			// split about median
	annMedianSplit(pa, pidx, n, cut_dim, cut_val, n_lo);
    }
}
